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Abstract. We discuss the transition from quasi-circular inspiral to plunge of a system of 
two nonrotating black holes of masses mi and ni2 in the extreme mass ratio limit mi?rt2 -C 
(mi +m2) 2 . In the spirit of the Effective One Body (EOB) approach to the general relativistic 
dynamics of binary systems, the dynamics of the two black hole system is represented in terms 
of an effective particle of mass /i = m\m2 / (mi + 7712) moving in a (quasi-)Schwarzschild 
background of mass M = mi +m2 and submitted to an C(/i) radiation reaction force defined 
by Pade resumming high-order Post-Newtonian results. We then complete this approach by 
numerically computing, a la Regge-Wheeler-Zerilli, the gravitational radiation emitted by such 
a particle. Several tests of the numerical procedure are presented. We focus on gravitational 
waveforms and the related energy and angular momentum losses. We view this work as a 
contribution to the matching between analytical and numerical methods within an EOB-type 
framework. 
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1. Introduction 

The last months have witnessed a decisive advance in Numerical Relativity, with different 
groups being able to simulate, for the first time and by using different techniques, the merger 
of two black holes of comparable masses (without or with initial spin) [1, 2, 3, 4]. Since 
such binary black holes systems (of a total mass^ 3OM ) are believed to be among the 
most promising sources of gravitational waves for the ground based detectors like LIGO and 
VIRGO, this breakthrough raises the hope to have, for the first time, a reliable estimate of 
the complete waveform by joining together Post-Newtonian (PN) and Numerical Relativity 
results. We recall that PN techniques have provided us with high-order results for describing 
the motion [6, 7] and radiation [8, 9] of binary systems, and that further techniques have been 
proposed for resumming the PN results [10, 11, 12], thereby allowing an analytical description 
of the gravitational waveform emitted during the transition between inspiral and plunge, and 
even during the subsequent merger and ringdown phases. We now face the important task of 
constructing accurate complete waveforms by matching together the information contained in 
Post-Newtonian and Numerical Relativity results. We view the present work as a contribution 
towards this goal (for a recent first cut at this problem see [13, 14, 15]). 

The work we present here belongs to a scientific lineage which started with the 
pioneering works of Regge and Wheeler [16], Zerilli [17], Davis, Ruffini, Press and Price [18] 



Binary black hole merger in the extreme mass ratio limit 



2 



and Davis, Ruffini and Tiomno [19]. References [18, 19] studied the gravitational wave 
emission due to the radial plunge (from infinity) of a particle onto a Schwarzschild black 
hole. This was thought of as a model for the head-on collision of two black holes in the 
extreme mass ratio limit. The gravitational wave energy spectrum [18] and waveform [19] 
were obtained. Davis, Ruffini and Tiomno [19] pointed out that the first part of the waveform 
could be described in terms essentially of a flat space "quadrupole formula" for a test particle 
following a Schwarzschild dynamics (Ruffini-Wheeler approximation [20]), while the last 
part of the waveform was dominated by exponentially damped harmonic oscillations. These 
damped oscillations were interpreted by Press [21] (see also Vishveshwara [22]) as vibrational 
modes (now called quasi-normal modes) of a Schwarzschild black hole. The perturbative 
formalism employed in these works has been later shown to be expressible in a gauge 
invariant manner [23, 24, 25, 26]. The case of a particle plunging from infinity with nonzero 
angular momentum has been discussed by Detweiler and Szedenits [27] and Oohara and 
Nakamura [28] by means of the curvature perturbation formalism of Teukolsky [29]. The 
radial plunge problem has been recently extended to a plunge from a finite distance [30, 31], 
with particular emphasis on the effect of the choice of initial data. 

However, none of the above works has studied the transition from the quasi-circular 
adiabatic inspiral phase to the plunge phase in extreme-mass-ratio binary black hole systems. 
The reason is that the original Regge-Wheeler-Zerilli test-particle approach seems to require 
the test particle to follow an exact geodesic of the Schwarzschild background. Here we shall 
bypass this stumbling block by appealing to some results of PN theory. In particular, the 
Effective One Body (EOB) approach to the general relativistic two body dynamics, which has 
been recently proposed to study the transition from inspiral to plunge in the comparable- 
mass case [11, 12, 34, 35, 36, 37], describes the dynamics of a binary system in terms 
of two separate ingredients: (i) a Hamiltonian Heob(M,/i) describing the conservative 
part of the relative dynamics, and (ii) a non-Hamiltonian supplementary force !Feob(M, (i) 
approximately describing the reaction to the loss of energy and angular momentum along 
quasi-circular orbits. Following a prescription suggested in [10] the badly convergent [38] PN- 
Taylor series giving the angular momentum flux is resummed by means of Pade approximants. 

Our general motivation for studying this problem is: (i) to gain information on black 
hole plunges in a regime (p <C M) that is not yet accessible to full numerical simulations, 
and (ii) to learn how to match analytical results to numerical ones in a situation where it 
is relatively easy to perform many, controllable high-accuracy numerical simulations. This 
paper will be mainly devoted to the discussion of our numerical framework; see [40] for a 
thorough comparison between analytical and numerical results. 



2. Relative dynamics of extreme-mass-ratio binary black holes 

In the EOB approach the relative dynamics of a binary system of masses mi and m,2 is 
described by a Hamiltonian Heob(M 7 fi) and a radiation reaction force Teob{M, /t), where 
M = mi + m 2 and \i = mim2/M. In the general comparable-mass case Heob has the 

structure Heob(M, ^) = M\J 1 + 2v(H u — 1) where v = [ijM = mim2/(mi + TO2) 2 
is the symmetric mass ratio. In the test mass limit that we are considering, v <C 1, we can 
expand Heob in powers of v. After subtracting inessential constants we get a Hamiltonian 
per unit (p) mass H = \im v ^ (H — const.)//! = lim^o H v of the form 

* = /«"( 1 + 4> + *)- <n 
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Here we have introduced the dimensionless variables f = r/M, p r = P r = P r /[i and 
p v = P v /M = P v /(/j,M). The functions A(f, v) and B(f, v) entering the Hamiltonian H v 
are metric coefficients entering the effective one body metric 

As 2 = g^dx^dx" = -A(f, v)At 2 +B{f, v)Ar 2 +r 2 (A6 2 + sin 2 6Aip 2 ) .(2) 
In the limit v — > in which we shall use (1) the effective metric functions A(f, v) and B(f, v) 
reduce to the well-known Schwarzschild expressions A(f,0) = (B(f, 0)) _1 = 1 — 2/f. 
Since the radiation generation process will be studied in terms of the Regge-Wheeler tortoise 
coordinate r* = r + 2M log(r/(2M) — 1), it is useful to have the Hamiltonian explicitly 
written in terms of r* and its conjugate momentum p Tt . This amounts to performing a 
canonical transformation (r,p r ) — > {r*,p r *)- The invariance of the action, p r ,Ar* — p r Af 
yields p r — (df * /df)p r . = A~ 1 p Vt , so that the new Hamiltonian function becomes 



H(f*,Pr.) 
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On the r.h.s. of the last equation we have introduced the second ingredient of the EOB 
approach: the radiation reaction force T v . The function /dis(w^) is the "factored flux 
function" of [10] scaled to the Newtonian (quadrupole) flux. Reference [10] has shown 
that the sequence of near-diagonal Pade approximants of /dis( w ^) exhibits, during the 
quasi-circular inspiral phase (f > 6), a good convergence toward the exact result known 
numerically [38] in the v — > limit. Here we shall use the v —> limit of the 2.5 PN accurate 
/Dis(wf) (i.e., equations (3.28)-(3.36) of [12] in the v — limit). Following the suggestion 
of [37], we expressed, in (8), T v as a function of the two variables u) and f with the hope 
that this expression remains approximately valid during the plunge (f < 6). We shall give 
evidence below that this hope is fulfilled. Note also that the flux used in constructing 
in (8), refers to the momentum flux at infinity. Reference [39] has shown that the flux down 
the event horizon (of the background Schwarzschild black hole of mass M) would correspond 
to a (here negligible) 4 PN effect. 



2.1. Quasi-circular orbits and transition from inspiral to plunge 

Following Section IV (A) of [12] and taking the v — > limit in all terms except the crucial 
T v = 0(v), we start the orbit at some given initial radius f and initial phase ip with 
corresponding momenta: 
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Figure 1. Left panel: Plunge relative orbit from r = 7M. The circle represented with a thick 
dashed line is the LSO at r = 6M. The thick solid line circle represents the light ring at 
r = 3M. Right panel: Radial velocity v r = f and orbital angular frequency ui = if versus 
coordinate time. 



where we explicitly have (with H = H (p r = 0)) 

(Ho) 2 V(r-3) 
B r = - 2 ~^A 2 (ff-^^, (13) 

Figure 1 shows the trajectory determined by this kind of initial f = 7, f = and v = 0.01 
(for the sake of convenience we fix M — 1 in the numerical calculations); after ~ 5 orbits the 
trajectory crosses the Last Stable Orbit (LSO) at f = 6 and then "plunges" (in a quasi-circular 
way u r <l, see inset in the right panel) into the black hole. It has to be noted that the orbital 
angular frequency Moo reaches a sharp maximum nearly at the light ring crossing (f = 3). 



3. Gravitational wave generation in the extreme mass ratio limit 

We now complete the description of an extreme-mass-ratio binary system by computing the 
gravitational wave emission driven by the system during the inspiral and the plunge. 

Let us recall that non-spherical linear metric perturbations around a Schwarzschild black 
hole can be expanded in scalar, vector and tensor spherical harmonics that are called even- 
parity [or odd-parity] if they transform under parity as (— I) 1 [or (— 1) £+1 ] and are decoupled 
because the black hole is non rotating. The seven even-parity metric multipoles [or the three 
odd-parity multipoles] can be combined together so that Einstein equations yield a wave-like 
equation for the even-parity [or odd-parity] gauge-invariant master variable [or "J^] 
with corresponding (Zerilli-Moncrief/Regge- Wheeler) potential V^°\ In Schwarzschild 
coordinates and in the presence of a general matter source driving the perturbation of the 
spacetime, they read 

«( 0) - a?.*& /0) + vt ,o) ^ o) - & . d4) 
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The well-known expressions of the potentials and the less well-known expressions of the 
sources for a general stress-energy tensor can be found in [25, 26]. We recall that from 
"J/^ and "J/^ one can easily obtain, when considering the limit r — > oo, the h + and h x 
gravitational-wave polarization amplitude, and the emitted power and angular momentum 
flux according to 



h + — \h y 



E 



1(1 + 


2)! 


V- 


2)! 


(*4 


-2)! 



-,(e) 
-Cm 



rim 



^ (l-2)\ 



— E 



(1 + 2)! 
''(*-2)! 



+ C.C. 



(15) 
(16) 
(17) 



where the overdot stands for time-derivative and _ 2 Y = -2Y (6, ip) are the s = 2 
spin-weighted spherical harmonics [41]. In our physical setting (y — > 0) the sourcing stress- 
energy tensor is, to leading order, that of a test particle following the (relative) dynamics 
described above. We can then use the gauge-invariant multipoles of the stress-energy tensor 
of a point particle as explicitly given in Appendix A of [26]. Since the Zerilli-Moncrief 
and Regge-Wheeler equations are solved on a r* grid, it is natural to express the source in 
terms of <5(r* — R*(t)) where R*(t) denotes the "particle tortoise coordinate", by contrast 
to a general "field coordinate" r*. In addition, since the particle dynamics is written using 
canonical variables, it is also natural to express the source in terms of the particle coordinates 
i?* and $ and their (rescaled) conjugate momenta P Tr and P v . This has the advantage [over 
the use of the r coordinate and 5(r — R(i)) ] of allowing one to push the evolution in time as 
much as one as one wants when r — ► 2M, without introducing spurious boundary effects that 
may spoil the waveforms (and in particular the late-time tail). The coordinate transformation 
r — > r* (with dr/dr* = 1 — 2M/r = A(r) ) implies the relations 

6(r-R(t)) =A(r)- 1 S{r*-R*{t)), (18) 

d r 5(r - R(t)) = - ™^5(r. - R.(t)) + J^^Jin - R*(t)) . (19) 

Straightforward algebra permits then to derive from the results of [26] the following explicit 
expressions for the source terms in (14) 
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+ (l - 5p) P v d r J(n - R*(t))j , (21) 

where A = £(£ + 1). In the above equations we have chosen to write the sources in the 
functional form 

S£i 0) = G { fl°\r,t)8{r* R*{t)) + F^° ] (r,t)d r J(r» R*(t)) , (22) 

with r-dependent (rather than i?(t)-dependent) coefficients G(r), F(r^. Note that the time 
dependence of F(r,t) and G(r,t) comes from the dependence on H, Pji t and P v . By 
exploiting the properties of the <5-function, we can also rewrite the sources in the functional 
form 

S%i 0) = G ( °l o) (RAt))5(r* - R*(t)) + F^°\R4t))d r J(n - R*(t)) , (23) 
where 



(24) 

Expressions (22) and (23) are mathematically equivalent, for a distributional source, but 
become numerically (slightly) different when the ^-function is approximated on the r* 
numerical domain by means of a narrow Gaussian. 

4. Numerical framework, tests and comparison with the literature 

4.1. Numerical procedure 

We solve the equations given in section 2 for the particle dynamics using a standard fourth- 
order Runge-Kutta algorithm with adaptive stepsize [42]. Then we insert the resulting 
position and momenta in the source terms <S^°^ using a Gaussian-function representation of 
S(r* — i?„ (t)) (see below). The corresponding Zerilli-Moncrief and Regge-Wheeler equations 
are then numerically solved, using standard techniques [42, 43], in the time domain (for each 
multipole (£, m)) by discretizing the axis ( € [r™ m , r™ ax ] ) with a uniform grid spacing 
Ar* . In particular, the solution can be computed either by means of a centred second-order 
finite differencing algorithm, or by the following implementation of the Lax-Wendroff method 
(that is the one usually preferred). In this second case, it is convenient to rewrite the wave- 
equations as a first-order flux conservative system (with sources) in the form d t U + d rr F = S, 
where we have defined the vector of "conserved quantities" U and the "fluxes" F as 

U = I ^e™ j F = f Sm J (25) 
where w = d t ^'fj°" 1 + d r _ t and the vector S is 

S = ( T/ (c/o). T .(o/o) q (e/o) ) , (26) 

Denoting by j the spatial grid-point index, n the time level and At the time step, the explicit 
time-advancing algorithm reads 

At At 2 
U ? +1 = U ? - 2A^ t F ? +1 ~ F ^ + 2A7f ^ ~ 2F ^ + ^ A + A * S " ' (2?) 

where we introduced the matrix A = diag(l, —1). This algorithm is stable under the standard 
Courant-Friedrichs-Levy [43] condition Ai/Ar* < 1. We typically use Ai/Ar» = 0.9. 
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On the boundaries we impose ingoing (at r™ m ) and outgoing (at r™ ax ) standard Sommerfeld 
conditions. This efficiently suppress reflections from the boundaries apart from tiny effects 
that can modify (if we do not use sufficiently large grids) the power-law tail at the end of 
the black hole ringdown phase. [Improved non-reflecting boundary conditions for Zerilli- 
Moncrief and Regge-Wheeler equations have been recently discussed [44]]. 



4.2. Approximating the 5-function 

We approximate the (5-function that appear in the source terms by a smooth function 5 a (r* ) 
(and d rr 5(r*) by 9 r „(5 CT (r*)). We usej 



6(n (*))-> -4=exp 

<7V27T 



(r, -R*(t)) 



2 1 



2a 2 



(28) 



with a > Ar*. In practice a ~ Ar* works well thanks to the effective averaging entailed 
by the fact that (t) is not restricted to the r* grid, but varies nearly continuously on the 
r* axis. In the continuum limit (a ~ Ar» — » 0) the derivatives of both and 
would be discontinuous at the location of the particle r* = generating numerical noise 
(Gibbs phenomenon) if standard numerical methods are used. However, the smoothing of 
<5(r* — 7?* (t)), together with the use of a numerical method like Lax-Wendroff (that has a bit 
of numerical dissipation built in) avoids any problems related to high frequency oscillations 
and provides us with a clean and stable evolution §. The important scale in our problem 
is M which determines the "width" of the Zerilli potential. We found that a resolution of 
a ~ Ar* = 0.01M is small enough to ensure numerical convergence to the continuum 
solution. 



4.3. Numerical tests: circular orbits and radial plunge 

We have tested the reliability of our code by recomputing results for simple geodesic 
trajectories {T v = 0) that have been already obtained in the literature using both time-domain 
and frequency domain approaches, as well as different treatments for the particle and different 
expressions for the sources. 

First of all we discuss the case of circular orbits. Table 1 lists the energy and angular 
momentum fluxes (at r Q b s = 1000M) up to the t = 4 multipole for an orbit of radius 
r = 7.9456. This permits a direct comparison with the results of Martel [32] (see also [33]), 
that we include in the table as well for the sake of completeness. For this test run we consider a 
(relatively) coarse resolution of Ar* = 0.02M; this is enough to have a very small difference 
with Martel results for the 1 — 2 multipoles, but the agreement obviously worsens for higher 
multipoles, for which it is necessary to increase the resolution. That is the reason why, when 
discussing the real plunge phenomenon driven by radiation reaction we shall use resolutions 
up to Ar* = 5 x 10~ 3 to properly capture the behaviour of the highest multipoles. 

Our second test concerns the waveform of a particle plunging radially into the black hole 
from a finite distance (where it starts at rest). This problem has received some attention in the 
recent literature [30, 31], thereby extending (by means of both frequency-domain and time 
domain approaches) the pioneering analysis of Davis et al. [18, 19] of the early 70s for a 
particle plunging from an infinite distance. 

\ Note that this smoothing of the distributional stress-energy tensor into a formally extended one only concerns the 
calculation of the waveform. The dynamics, equations (4)-(8), is that of a (5-function stress-energy tensor. We use 
this smoothing only as a numerical technique, and we have checked that we were in the convergent regime. 
§ Note that the use of a method without dissipation, like a standard leapfrog, introduces high frequency noise if a is 
too small. For this reason we prefer to use the Lax-Wendroff technique for most of the simulations. 
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Table 1. Energy and angular momentum fluxes extracted at r b s = 1000M for a particle 
orbiting the black hole on a circular orbit of radius r = 7.9456. Comparison with the results 
of Mattel [32]. 
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Figure 2. Test of the code: waveform (on a logarithmic scale in the inset) for a particle 
plunging radially on the black hole along the z-axis from r = 10M 



We selected the same kind of initial data as described in references [30, 31], so as to have 
an initial profile of \I>^ that is conformally flat. However, since we have slightly "smeared" 
the (5-function, we cannot use the "discontinuous" analytical initial data of references [30, 31]. 
We numerically solve the Hamiltonian constraint by writing it as a tridiagonal system that is 
then solved using a standard LU decomposition. Figure 2 displays (for the case I = 2) the 
waveform generated by a particle plunging into the black hole along the z-axis, starting from 
rest at r = 10M. It has been extracted at r b s = 500M and is shown versus the retarded 
time u = t — r° bs . The master function ty^o nas been multiplied by a factor two (linked to a 
different choice of normalization in [30, 31] ) to facilitate the (very satisfactory) comparison 
with the top-right panel of figure 4 in [31] or the top-left panel of figure 6 in [30]. We used 
a resolution of Ar* = 0.01 M with r* e [-800M, 1800M]. This avoids the influence of 
the boundaries and allows one to capture the late-time non-oscillatory decay after the QNM 
ringdown phase (see the inset in figure 2). 
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Figure 3. Dominant £ = m = 2 even-parity Left panel and and I = 2, m = 1 odd-parity 
Right panel gravitational wave multipoles generated by the quasi-circular plunge of a particle 
with v = 0.01 initially at r = 7M. 



5. Transition from inspiral to plunge: waveforms 

Now that we have shown that our code can reproduce with good accuracy known results for 
geodesic orbits, let us break new ground by considering the waveform emitted during the 
transition from inspiral to plunge. We refer here to the dynamics described in section 2.1; 
i.e., the particle with v = 0.01 is initially at r = 7A1 with dynamical initial data defined by 
the quasi-circular (adiabatic) approximation described there. We also need initial data for the 
field perturbations: ^ < fj l °\ d^f^ . To do so, one should in principle select the solution 
of the perturbed Hamiltonian and momentum constraints which corresponds to the physically 
desired "no-incoming radiation condition". However, we took a more pragmatic approach. 
As it is often done in the literature, we impose X P^°' > = dt^^J^ = at t = 0. This "bad" 
choice of initial data produces an initial burst of unphysical radiation, but after a while the 
gravitational waveform is driven by the motion of the source so that it is sufficient not to 
include the early part of the waveform in the final analysis. Figure 3 depicts the dominant 
1 = 2 waveforms {m = 2 and m = 1 multipoles) extracted at r b s = 250A/. For the 
sake of comparison we normalize the waveforms to the mass /i (see below a discussion of 
the approximate universality of this scaled waveform). The numerical grid we adopted for 
this simulation was r„ £ [— 1200M, 1200M] with resolution Ar* = 0.01AZ. As we were 
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Table 2. Total energy and angular momentum extracted at r t, s = 250M for the "plunge 
phase" (r < 5.9865M) of a dynamics with v = 0.01. 
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mentioning above, after an initial unphysical burst of radiation the gravitational waveforms 
are driven by the motion of the particle during the quasi-circular inspiral and plunge. From 
the inspection of the particle dynamics we described above, we know it crosses the LSO at 
u/(2M) ~ 240 and passes through the light ring (where the gravitational perturbation driven 
by the source is filtered by the peak of the potential) at u/(2M) ~ 300 (where both the 
orbital frequency u and the radial velocity v r have a maximum). Both and \I>^ show a 
progressive increase in frequency and amplitude until uj (2M) ~ 300, where a maximum of 
amplitude is reached. Then follows a QNM dominated ringdown phase. 

To validate our results, we study the convergence of the waveforms by considering 
resolutions Ar» = (0.02, 0.01, 0.005). For the dominant I = m = 2 waveforms we found a 
convergence rate (3 ~ 1.6, defined from oc Ar* . Here A^P is defined as the discretized 
root mean square between the waveform at the resolution Ar* and the highest resolution one. 
Generally speaking, a resolution Ar* = 0.01 is sufficient to determine the energy and angular 
momentum radiated with an accuracy that is of the order of 1% (or better) for the quadrupole 
modes. Having the same accuracy for higher modes needs to increase the resolution; however, 
we have verified that this is sufficient for having an accuracy (for the energy) that is not worse 
than 6% for the other multipoles up to I = 4 (the most sensitive one being t = 4, m = 0). 

To give some numbers we computed the total energy and angular momentum loss during 
what might (approximately) be called the plunge phase (that is, after the crossing of the LSO 
at 6M). More precisely, we selected the part of the waveform for u/(2M) > 240 which 
corresponds to radii f < 5.9865. Table 2 lists the partial multipolar energy and angular 
momentum losses up to I = 4. These values were computed with resolution Ar* = 0.01 so 
that the accuracies are of the order of 1% for t = 2, 2 - 4% for I = 3 and of 5 - 6% for t = 4. 
When one sums all the multipoles, the total energy emitted is found to be MA£//i 2 ~ 0.51 
and AJ//i 2 ~ 4.3. Only for the sake of comparison, let us mention that this "plunge" 
MAE/^i 2 is about 50 times larger than the "radial" (Davis et al. [18, 19]) one (summed 
up to I = 4) which amounts to MAE/^i 2 = 0.0104 (for a plunge from infinity). 

An important consistency check of our approach ( which assumes the specific PN- 
resummed radiation reaction force (8) ) is to compare the angular momentum loss assumed 
in the dynamics, i.e — (dP v / dt) / (/J.M) = —T^, to the angular momentum flux radiated 
by the multipolar gravitational waves at infinity (up to I = 4). The left panel of 

figure 4 displays — (dP v / dt) / (fiM) (dashed line) versus (dJ/dt) gv/ /(fiM) (solid line) up to 




Figure 4. Left panel: angular momentum flux computed from the gravitational waveform 
compared to the angular momentum loss assumed in the dynamics with v = 0.01. Right panel: 
Instantaneous gravitational wave frequency for two values of v. The signal is " approximately 
universal" after u/(2M) ~ 280, which roughly corresponds to r < 5.15 ("quasi-geodesic 
plunge"). After u/(2M) ~ 300 starts the QNM phase. 



roughly the crossing of the light ring. The very good agreement between the two curves is 
a confirmation of the good convergence of the Pade-resummed radiation reaction force F v , 
equation (8), toward the exact result. This confirmation goes beyond the tests of [10] which 
were limited to the inspiral phase f > 6. The retarded times 240 < u/(2M) < 290 in figure 4 
correspond to the plunge phase [the merger (matching) occurring at u/(2M) ~ 300]. 

The question that one can ask at this point is how much the numbers of table 2 
are universal. Let us first recall that [12] has shown that the transition between inspiral 
and plunge was taking place in a radial domain around the LSO which scaled with v as 
r - 6M ~ ±1.89MV 2 / 5 , with radial velocity scaling as v r ~ -0.072^ 3/5 . The specific 
power 2/5-th appearing in the radial scaling yields, in the numerical case v = 0.01 that we 
consider in most of this paper, r— 6M ~ ±0.30M. This means that, formally, when v = 0.01, 
it is only when 6M—r ^S> 0.3M that we can expect the plunge dynamics to become universal, 
in the sense that it will approach the geodesic which started from the LSO in the infinite past 
with zero radial velocity. If we wanted to reach this universal behaviour just below the LSO 
(say, for r < 5.98) we would need to use much smaller values of v (say v < 10~ 5 ). For 
larger values of v we can only hope to see an approximate convergence to universal behaviour 
near the end of the plunge. This is illustrated in the right panel of Figure 4. This figure plots, 
versus retarded time u, the instantaneous gravitational wave frequency Mw^, where (3 
imaginary part) 

M <™ - ^ Ufc (29) 

for v = 0.01 (solid line) and v — 0.001 (dashed line). At early times, that is when the particle 
moves along a quasi-circular orbit, one has Muj 1 ™ ~ mMu) — 2Mlq\ i.e., the double of the 
orbital frequency. Then one has the transition regime mentioned above around the crossing 
of the LSO, i.e., for u/(2M) ~ 240. The quasi-universal, quasi-geodesic plunge mentioned 
above starts afterwards around uj (2M) ~ 280, i.e., r ~ 5.15M (for v = 0.01) 
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Figure 5. Consistency checks of the sources. Left panel: relative difference between 
waveforms obtained with different expressions for S22 ■ Right panel: check of the smallness 
of the effects linked to the terms oc P v and oc H in the odd-parity source Soi from (21). 



Then, after crossing the light ring, one enters another universal phase: the QNM one. In 
this QNM phase the gravitational wave angular frequency saturates and oscillates around the 
value Mui^j = 0.3736715 (highlighted in the inset). The fine structure of this gravitational 
wave frequency plot will be further discussed in [40]. 

We conclude this section by briefly quoting two consistency checks of the matter source 
that we have performed. The first check concerns the consistency between the two different 
functional forms (22) or (23) that the sources can take. As said above these two forms would 
be mathematically equivalent for a real (5-function, but will differ when using the Gaussian 
approximation (28). The left panel of figure 5 shows the relative difference between two 
gravitational wave moduli 1^22 I (dotted line): one of them computed using (22) and the other 
using (23) with Ar* = 0.02 and a — Ar*. This left panel also shows the difference between 
the gravitational wave frequencies Mw^ (solid line) computed with the two different sources. 
The relative difference ~ 10~ 5 on both makes us confident that the "smeared" (5-function is a 
very good approximation to the actual ^-function. [Let us also comment, in passing, that we 
did a convergence check (for £ — m = 2) by comparing waveforms obtained by means of 
Gaussians (with a = Ar* and a — 2Ar*), finding a relative difference of the same order of 
magnitude ~ 10~ 5 ]. 

Finally, we performed a consistency check regarding the explicit time derivative 
appearing in the coefficient of S(r* — i?*(t)) in the odd-parity source S$ from (21). When 
expanded by Leibniz rule, this time derivative generates three terms. The term proportional 

to Pr, would be present along a geodesic motion. By contrast, the terms proportional to P v 

and H are proportional to v and therefore vanish in the geodesic motion limit. These last two 
terms generate a correction of order 0(n)0(v) oc v 2 in the source. Such terms are formally 
negligible in the extreme mass ratio limit that we consider here. In the right panel of figure 5 
we have checked to what extent these terms are numerically small. 
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Figure 6. Comparison between the gravitational-wave phase for I = m = 2 obtained from 
the Zerilli-Moncrief equation or from analytically matching a 3PN improved quadrupole-type 
formula (factor F22 Pade) to a superposition of quasi-normal modes. 



6. Conclusions 

As we briefly mentioned in the Introduction, beyond solving for the first time, within a 
certain approximation, the problem of the plunge in the extreme mass ratio limit, the most 
important aim of this work was to build "actual" numerical waveforms to be then compared 
with analytical ones, within the EOB framework and philosophy. 

We recall that the basic idea of the EOB framework is to produce quasi-analytical 
waveforms by patching together a quadrupole-type waveform during the inspiral and plunge 
to a QNM-type waveform after merger. In [40] we will use the numerical tools presented here 
to measure to what extent these EOB-type waveforms can approximate actual waveforms 
of binary black hole merger in the extreme mass ratio limit. We give an example of 
this numerical-analytical comparison in figure 6. This figure compares the phases of 
the gravitational waveforms obtained by two different methods, one numerical and the 
other semi-analytical. In the first method, the phase has been computed by integrating 
in time the instantaneous gravitational wave frequency given by (29), with obtained 
by numerically solving the Zerilli-Moncrief equation. The second method computes an 
approximate waveform in the following way: Before crossing the light ring one uses a (Pade- 
resummed) 3PN-improved quadrupole-type formula to compute the waveform from the EOB 
dynamics. After crossing the light ring, the previous quadrupole-type signal, taken in the 
quasi-circular (QC) approximation, is matched to a superposition of the first five QNMs of the 
black hole. Then one computes the phase of this matched analytical waveform by integrating 
the corresponding instantaneous gravitational wave frequency. The agreement between the 
"actual" phase and the "effective one body" phase is impressively good: it turns out that the 
maximum difference between the two is less than 3% of a cycle. 
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